## Linear Programming upper bound

from PyM import *

def LP(n, d, q = 2, R = []):
    if q==2: return LPb(n,d,R)  # special processing when q=2
    K = [krawtchouk(k,n,q) for k in range(1,n + 1)]
    A = [[evaluate(k,0)] + [evaluate(k,i) for i in range(d,n+1)] for k in K]
    A = -matrix(A)
    alpha = symbol('a',n)
    alpha = [1] + alpha[d - 1:n]
    alpha = vector(alpha,'c')
    RT = A*alpha
    RT = [r for r in RT] + R
    M = sum(alpha)
    m1,m2 = simplex(M,RT)
    a0 = symbol('a0')
    return floor(m1),{**{a0:1},**m2}

def LPb(n,d,R=[]):
    if odd(d): n = n + 1; d = d + 1
    K = [krawtchouk(k,n) for k in range(1,n//2 + 1)]  
    A = [[evaluate(k,0)] + [evaluate(k,i) for i in range(d,n+1,2)] for k in K]
    A = -matrix(A)
    alpha = symbol('a',n)  
    alpha = [1] + alpha[d - 1:n:2]
    alpha = vector(alpha,'c')
    RT = A*alpha
    RT = [r for r in RT] + R
    M = sum(alpha)
    m1,m2 = simplex(M,RT)
    b = floor(m1)
    if odd(b):
        beta = stack([1-1/(b>>Q_)],alpha[1:])
        RT = A*beta
        RT = [r for r in RT] + R
        m1,m2 = simplex(M,RT)
    a0 = symbol('a0')
    return floor(m1),{**{a0:1},**m2}

[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14] = symbol('a',14)
show(LP(13,5))
show(LP(14,6))
[a10,a12] = symbol('a10','a12')      
A = LP(13,6, R=[a10+4*a12-4])   
show(A)
[a6,a8] = symbol('a6','a8')
B = LP(9,4,R=[a8-1,a6+4*a8-12])
show(B)  

# Golay [11,6,5]_3

show(LP(11,5,3))

show(LP(12,6,3))

# Golay [23,12,7]: LP(23,7) = LP(24,8)

show(LP(24,8))